%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;

% ESTIMATING Full MODEL PARAMETERS			       		;
% T and F0 are now known at time=0                      ;
% T and F0 can now be correlated                        :
% Multi-Stage Optimization                              ;
% Oct 20, 2014 -- Paulina Oliva & Chris Severen         ;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%;

clear all;
warning off all;
format compact
format shortg

load resultsFMIN.mat

theta_hat = theta_hat2;

%*******************************************************;
% Execution and Optimization Options                    ;
%*******************************************************;


% rp.parl    = 1;     % Run in parallel: 0 - serial, 1 - parallel
% rp.cores   = 12;    % Only used if parallel==1
% rp.k       = 1500;
% rp.m       = 100;
% rp.lam     = 0.5;     % Smoothing parameter
% rp.tol     = 1e-15; % Tolerance in LL
% rp.theseed = 25000; % RNG initilization
rp.cc      = Inf; % Convergence criterion

params     = [];
stage      = [];

if rp.parl==0
    % Note: If run serial, it is faster to generate everything once, as
    % done here. If run in parallel, it is faster to generate these locally
    % than using up transfer bandwidth. This creates everything and
    % collects it into the structure snPass.
    
    rng(rp.theseed);
    snPass.s       = rand(obs*(3*rp.k + rp.m),1);
    
    % These two lines create a (k x m x 51) matrix, where the first floor is a
    % (k x m) matrix of ones, the second a (k x m) matrix of twos, up to 51;
    % Nv is just N repeated k x 51 times;
    % the fourth and fifth are true false matrices for N vs. Nbar;
    Ntr            = permute(N,[3,1,2]);
    snPass.N3D     = Ntr(ones(rp.k,1),ones(rp.m,1),:);
    snPass.Nv      = N(ones(rp.k,1),:);
    snPass.Ngeq3D  = (snPass.N3D>=Nbar);
    snPass.N03D    = (snPass.N3D>0);
    snPass.Ngeqv   = (snPass.Nv>=Nbar);
    snPass.N0v     = (snPass.Nv>0);
elseif rp.parl==1
    snPass.s       = 0;
end

%*******************************************************;
% Setting up values                                     ;
%*******************************************************;

hv = [0.0001,0.0005, 0.001:0.002:0.01, 0.01:0.02:0.1, 0.1:0.2:1, 2,3,4,5];


gprime(1) = 1;
gprime(2) = exp(theta_hat(2));
gprime(3) = 2*exp(-theta_hat(3))/((exp(-theta_hat(3))+1)^2);
gprime(4) = exp(theta_hat(4));
gprime(5) = exp(theta_hat(5));
gprime(6) = 1;
gprime(7) = 1;
gprime(8) = 1;
gprime(9) = 1;

gprime    = diag(gprime);

if rp.parl==1
    matlabpool('open',rp.cores)
end

for k=1:length(hv)
    
    thetamat = theta_hat(ones(1,18),:);
    
    thetamat(1:2,1)   = linspace(theta_hat(1)-hv(k),theta_hat(1)+hv(k),2)';
    thetamat(3:4,2)   = linspace(theta_hat(2)-hv(k),theta_hat(2)+hv(k),2)';
    thetamat(5:6,3)   = linspace(theta_hat(3)-hv(k),theta_hat(3)+hv(k),2)';
    thetamat(7:8,4)   = linspace(theta_hat(4)-hv(k),theta_hat(4)+hv(k),2)';
    thetamat(9:10,5)  = linspace(theta_hat(5)-hv(k),theta_hat(5)+hv(k),2)';
    thetamat(11:12,6) = linspace(theta_hat(6)-hv(k),theta_hat(6)+hv(k),2)';
    thetamat(13:14,7) = linspace(theta_hat(7)-hv(k),theta_hat(7)+hv(k),2)';
    thetamat(15:16,8) = linspace(theta_hat(8)-hv(k),theta_hat(8)+hv(k),2)';
    thetamat(17:18,9) = linspace(theta_hat(9)-hv(k),theta_hat(9)+hv(k),2)';
    
    lik      = zeros(1,18);
    scorevec = zeros(obs,18);
    
    %*******************************************************;
    % Function call, start min                              ;
    %*******************************************************;
    
    for i=1:18
        [lik(i),scorevec(:,i)] = Full_Likelihood(thetamat(i,:),DP,I_N,A,R,TG,MON,params,rp,snPass,stage);
    end
    
    [grad,hess1,hess2] = Full_HessScore(thetamat,scorevec);
     
    varcovar1 = gprime * inv(hess1) * gprime' / obs;
    varcovar2 = gprime * inv(hess2) * gprime' / obs;
    sevec1    = sqrt(diag(varcovar1));
    sevec2    = sqrt(diag(varcovar2));
    %output2vec = [coefftran(theta_hat,0);sevec'];
    
    seresults(k).grad  = grad;
    seresults(k).hess1 = hess1;
    seresults(k).hess2 = hess2;
    seresults(k).hv    = hv(k);
end

if rp.parl==1
    matlabpool close
end

save resultsSE.mat
